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Abstract. We have implemented an Adaptive Mesh Refinement criterion explicitly designed to increase spatial resolution 
around discontinuities in the velocity field in ENZO cosmological simulations. With this technique, shocks and turbulent eddies 
developed during the hierarchical assembly of galaxy clusters are followed with unprecedented spatial resolution, even at large 
distances from the clusters center. By measuring the spectral properties of the gas velocity field, its time evolution and the 
properties of shocks for a reference galaxy cluster, we investigate the connection between accretion processes and the onset of 
chaotic motions in the simulated Inter Galactic Medium over a wide range of scales. 
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1. Introduction 

The intergalactic medium (IGM) in galaxy clusters is likely tur- 
bulent, at some level: this is claimed from several independent 
theoretical and numerical approaches (e.g. Bryan & Norman 
1998; Ricker & Sarazin 2001; Brunetti et al.2001; Inogamov 
& Sunyaev 2003; Dolag et al.2005; Subramanian et al.2006; 
Vazza et al.2006; Brunetti & Lazarian 2007; Nagai et al.2007; 
Iapichino & Niemeyer 2008). A number of observational ev- 
idences has also been published in the last few years. Using 
a mosaic of XMM-Newton observations of the Coma cluster, 
Schuecker et al.(2004) obtained spatially-resolved gas pressure 
maps which indicate the presence of a significant amount of 
turbulence, with a spectrum of the fluctuations consistent with 
a Kolmogorov turbulence. Additional evidences of turbulent 
motions inside nearby galaxy clusters came from the observa- 
tion of pseudo-pressure fluctuations in Abell 754 using XMM 
(Henry et al.2004) and from the non detection of resonant scat- 
tering in the Perseus cluster (Churazov et al.2004). Also studies 
of Faraday Rotation allow a complementary approach and sug- 
gest that the IGM magnetic field is turbulent on a broad range 
of scales (Murgia et al.2004; Govoni et al.2006; Ensslin & Vogt 
(2006). 

Detailed X-ray analysis performed in nearby cool-core 
galaxy clusters (e.g. Fabian et al.2003; Churazov et al.2004; 
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Graham et al.2006; Ota et al. 2006) suggest that the turbulent 
velocity field is subsonic at the scale of the cluster cores. Also, 
limits to the amount of turbulence in the IGM were recently de- 
rived by Churazov et al.(2008), suggesting that the amount of 
non-thermal pressure within ~ 50 kpc from the central galaxies 
in Perseus and Virgo clusters cannot exceed ~ 10-20 per cent 
of the thermal energy budget at the same radius. 

In addition, the phenomenology of diffuse radio halo emis- 
sion suggests a scenario in which turbulent MHD modes, ex- 
cited during cluster mergers, may re-accelerate the relativistic 
emitting particles (e.g. Ferrari et al.2008; Brunetti et al.2008; 
Cassano 2009 and references therein). Remarkably, the inter- 
play between Cosmic Rays (CR) and turbulent magnetic fields 
may drive still poorly explored plasma processes that may po- 
tentially affect our simplified view of the IGM (Subramanian 
et al. 2006, Schekochihin et al.2007, Brunetti & Lazarian 2007, 
Guo & Oh 2008). From the theoretical point of view, turbulence 
can be injected in the IGM by several mechanisms: plasma in- 
stabilities, cluster mergers and shock waves, wakes of galaxies 
moving into the IGM, outflows from AGNs hosted in the center 
of galaxy clusters and galactic winds. The total energy budget 
in form of turbulent motions inside galaxy clusters, as well as 
their distribution and their connection with cluster dynamics 
and non gravitational process in galaxy clusters are sill open 
fields and cosmological numerical simulations are potentially 
able to provide a great insight in the characterization of the 
above phenomena. 

Early Eulerian numerical simulations of merging clusters 
(e.g., Bryan & Norman 1998; Ricker & Sarazin 2001) provided 
the first reliable representations of the way in which turbulence 
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Table 1. Main characteristics of the runs. "D" stands for AMR 
triggered by gas/DM over-density, while "V" stands for AMR 
triggered by velocity jumps. A is the peak gas spatial resolution. 
6 specifies the value adopted to trigger AMR, see Sec. 3 for 
explanations. 



ID 


grid 


M dm [MJh] 


A [kpc] 


AMR 


v256-4 


256 3 


6.76 • 10 s 


18 


D+V(<5 > 10) 


v256-3 


256 3 


6.76 • 10 s 


36 


D+V(<5 > 3) 


vl28-3 


128 3 


5.39 ■ 10 9 


36 


D+V(<5 > 3) 


v64-3 


64 3 


4.32 • 10'° 


36 


D+V(<5 > 3) 


dl28 


128 3 


5.39 • 10 9 


36 


D 


V128-10 


128 3 


5.39 • 10 9 


36 


D+V(<5 > 10) 


vl28-l 


128 3 


5.39 • 10 9 


36 


D+V((5 > 1) 


vl28-z2 


128 3 


5.39 • 10 9 


36 


D+V(<5 > 3,z > 2) 



is injected into the IGM by merger events. More recently, high 
resolution Lagrangian (Dolag et al.2005, Vazza et al.2006) and 
Eulerian simulations (Nagai et al.2007, Iapichino & Niemeyer 
2008) found that a sizable amount of pressure support (i.e. ~ 
10 percent of the total pressure inside 0.5/?, 7 >) in the IGM is 
sustained by chaotic motions. Also, the amount of turbulent 
energy delivered by mergers and accretions is found to scale 
with the thermal energy of simulated galaxy clusters (Vazza et 
al.2006). 

Despite the tremendous capability that Lagrangian simula- 
tions have in resolving the smallest structures within galaxy 
clusters, they may suffer of serious limitations in modeling 
fluid instabilities, mostly because of the effects played by the 
artificial viscosity employed to solve hydro equations (e.g. 
Agertz et al.2007; Tasker et al.2008; Mitchell et al.2009). 
Therefore, the use of an Eulerian scheme free of artificial 
viscosity, as the Piecewise Parabolic Method adopted in the 
ENZO, can provide an important insight in all the above points. 
On the other hand Eulerian schemes with fixed grid resolution 
are typically limited by their low spatial resolution so that the 
application of Adaptive Mesh Refinement (AMR) techniques is 
mandatory to achieve adequate spatial detail in the simulations. 

In this work we present first results from the application 
of an Adaptive Mesh Refinement (AMR) to ENZO simula- 
tions, which allows to follow at the same time shocks and tur- 
bulent motions with unprecedented resolution up to large dis- 
tances from the galaxy cluster center. For the simulations pre- 
sented here, we assume a ACDM cosmology with parameters 
O = 1.0, Q BM = 0.0441, Q. DM = 0.2139, O a = 0.742, Hubble 
parameter h = 0.72 and a normalization of <x 8 = 0.8 for the 
primordial density power spectrum. 

2. Numerical Code and Setup 

ENZO is an adaptive mesh refinement (AMR) cosmological 
hybrid code highly optimized for supercomputing (Bryan & 
Norman 1997, 1998; Norman & Bryan 1999; Bryan, Abel, & 
Norman 2001, O'Shea et al.2004; Norman et al.2007). It cou- 
ples an N-body particle-mesh solver with an adaptive mesh 
method for ideal fluid-dynamics (Berger & Colella, 1989). 
ENZO adopts an Eulerian hydrodynamical solver based on the 
the Piecewise Parabolic Method (PPM, Woodward & Colella, 
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Fig. 1. Redshift evolution of the Dark matter plus gas mass 
{black), and of total gas mass {blue) inside the virial radius of 
the galaxy cluster studied in this work. 

1984), that is a higher order extension of Godunov's shock cap- 
turing method (Godunov 1959). The PPM algorithm belongs 
to a class of schemes in which an accurate representation of 
flow discontinuities is made possible by building into the nu- 
merical method the calculation of the propagation and interac- 
tion of non-linear waves. It is at least second-order accurate in 
space (up to the fourth-order, in the case of smooth flows and 
small time-steps) and second-order accurate in time. The PPM 
method describes shocks with high accuracy and has no need 
of artificial viscosity, leading to an optimal treatment of energy 
conversion processes, to the minimization of errors due to the 
finite size of the cells of the grid and to a spatial resolution close 
to the nominal one. In the cosmological framework, the basic 
PPM technique has been modified to include the gravitational 
interaction and the expansion of the Universe. 
We present here the simulation of a cubic volume of side 
75 Mpc starting from z = 30, and applying AMR within a sub- 
volume of side 1 .5Mpc, centered on a ~ 2 • 1O 14 M galaxy 
cluster. We re-simulate this volume under different configura- 
tions, as reported in Tab. 1 . The mass resolution of Dark Matter 
(DM) particles ranges from 6.76 • 10 8 M o (v256-3 and v256- 
4) to 4.32 ■ 10 10 M o (v64-3), corresponding to minimum root 
grid spatial resolutions from 292kpc to \.\12Mpc. The maxi- 
mum spatial resolution in the region where AMR is applied is 
A = 36kpc in all the simulations except for the case of v256-4, 
where A = \2>kpc. All runs are non radiative, and furthermore 
no treatment of reionization background due to AGN and or 
massive stars is modeled here. 

In all the above simulations, the galaxy cluster is formed 
through a major merger at 0.8 < z < 1, and visual inspection 
shows that its perturbed dynamical state stays till later epochs, 
due to further accretions. Fig[T] shows the redshift evolution 
of the total mass and of the gas mass inside the virial cluster 
region, measured according to the spherical over-density cen- 
tered on the density peak where the galaxy cluster forms. 
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Fig. 2. Gas density and temperature slices for the AMR re- 
gion of the vl28-3 run (upper panels), and of the dl28 run 
(lower panel). The gas density is normalized to the cosmo- 
logical critical gas density: rhocr = p cr = p cr ,ofb, where 
pcr,o = 3H(z) 2 /87tG is the cosmological critical density ( H(z) 
is the Hubble constant), and ft, is the cosmic baryon fraction. 
The side of the image is 7.5Mpc and the depth along the line 
of sight is 36kpc. The gas density is normalized to the critical 
density of the universe, rescaled by the cosmic baryon fraction. 

Computations described in this work were performed 
using the ENZO code developed by the Laboratory for 
Computational Astrophysics at the University of California in 
San Diego (http://lca.ucsd.edu). 



3. Adaptive Mesh Refinement technique for 
Turbulent Motions 

The first application of AMR to the study of turbulence in 
the inter stellar medium was reported in Kritsuk, Norman & 
Padoan (2006). Iapichino & Niemeyer (2008) applied a refine- 
ment criterion based on the gas velocity field (analyzing curl 
and divergence of velocity), in order to study turbulence in cos- 
mological ENZO simulations. Motivated by the above results, 
here we report on an exploratory study where a grid refine- 
ment scheme based on the analysis of one dimensional jumps 
in the velocity field is introduced in ENZO. In order to apply 
this method at full power to ENZO simulations we combine 
the implementation of the standard grid refinement criterion, 
customary adopted in cosmological simulations (based on over 
density), with a new grid refinement criterion based on the anal- 
ysis of the jump of velocity, Av, across cells. This choice en- 



sures that shocks can be studied with the highest available res- 
olution in simulated galaxy clusters, contrary to usual AMR 
runs, and at the same time this refinement scheme allows to 
increase the spatial resolution around turbulent features in the 
simulated galaxy clusters. 

In Teyssier (2002) 1-D tests are presented to assess the 
importance of refining shocks (according to a pressure crite- 
rion) in cosmological simulations. The major findings were 
that: a) relevant numerical instabilities do occur when simu- 
lated shocks move from a low resolution to a high resolution 
region (e.g. from a low density to a high density environment); 
b) since most of cosmological shocks move from high density 
(collapsing) regions to low density regions, explicitly refining 
on shocks can be safely avoided in the run time calculation of 
expanding accretion shocks (such as those developed in a stan- 
dard Zeldovich pancake collapse). However, in the high reso- 
lution simulation of galaxy cluster we present here, we expect 
significant departures from any idealized self-similar model of 
shocks (e.g. Molnar et al.2009), and the additional refinement 
on shock waves is a interesting option. 

In more detail, we propose to use the normalized 1-D ve- 
locity jump across 1-D patches in the simulation, 6 = |Av/v m | 
(where v m is the minimum velocity, in modulus, over the cells 
in the patch) to trigger grid refinement. Even if this method 
is highly simplified respect to that employed in Iapichino & 
Niemeyer (2008), in the next Sections we will show that it pro- 
duces a significant step forward, with no significant extra ex- 
pense of computational effort, in both the study of shock waves 
and the spectral characterization of the gas velocity field inside 
galaxy clusters. 

In particular we recursively analyze the velocity jumps 
across three adjacent cells at a given AMR level, and increase 
the resolution (by a factor 2 in cell size) for the cells of the patch 
whenever 6 is larger than a threshold value. At the same time, 
also the standard AMR method triggered by gas/DM over- 
density is applied (e.g. Norman et al.2007); the over-density 
threshold is set 6 P = 6p/p = 2 (where p can be either Dark 
Matter of gas density) for all runs. We notice that this threshold 
is smaller than what usually taken in similar works (e.g. 6 P — 4 
in O'Shea et al.2004; Nagai et al.2006; Iapichino & Niemeyer 
2008.) and thus typically much more volume is refined in the 
simulations presented here. 

We adopt as reference value 6-3 and allow for a number 
of AMR levels up to the maximum resolution of A = 36kpc. 
In one case, we also perform a run using the same setup of 
the v256-3 run, but allowing for one more AMR level (4 lev- 
els instead of 3), reaching the maximum resolution of \%kpc 
(v256-4). Finally, we present results for 6 = 10 (vl28-10) and 
6-1 (vl28-l), in order to assess the convergence of ourresults 
(Sec l43P31) . 

A reference simulation is also produced where only the gas/DM 
over-density criterion is used to trigger mesh refinements 
(dl28), along with a test run where the AMR criterion trig- 
gered by velocity jumps is added to the standard one only start- 
ing from z < 2 (vl28-z2). The latter run is designed in order 
to establish whether it is feasible to apply the novel mesh re- 
finement criterion starting only at later cosmic epochs, where 
clusters formation starts, saving some computational effort. 
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Fig. 3. Temperature maps for a central slice in the simulated 
AMR, at for different redshifts (z= 2, z= 1 and z=0.2) by using 
the standard AMR criterion (dl28 run, left panels), the new 
AMR criterion (vl28-3 center panels); the right panels show 
the cell by cell difference, as T new - T standard- The color table is 
as in Fig|2] right panels. 



4. Results 

4.1. Comparison with Standard AMR Runs. 

Fig|2] shows 2-D slices of gas density and temperature com- 
paring runs vl28-3 and dl28 at z = 0.1. Unlike in the standard 
mesh refinement triggered by gas/DM over-density, with the 
new AMR criterion shocks and chaotic motions are followed 
at the highest available resolution in the run, A = 36kpc, up to 
large (~ 3 — AMpc) distances from the cluster center. The dif- 
ference between the two approaches is remarkable at all stages 
in the evolution of the cluster, and especially in highlighting 
strong shock waves excited during the major merger event, as 
shown in the temperature maps of Figf3] 

In FigHJwe show the distribution of the volume occupied 
by cells at the different available AMR levels, comparing the 
results from vl28-3 and dl28 at z = 0.1. This shows that the 
application to the AMR criterion triggered by velocity jump 
does not cause any appreciable increase of memory expense in 
cosmological numerical simulations, compared to the adoption 
of the standard AMR criterion. Interestingly enough, although 
the volume occupied by cells at the highest AMR level is simi- 
lar in both runs (i.e. ~ 55 per cent of the AMR volume of side 
7.5 Mpc), we measure a slightly larger number of refined cells 
in the dl28 run (~ 1 — 2 per cent) at all AMR level. Differences 
in the distributions are consistent with the effect driven by the 
differences in the thermodynamics of the gas simulated with 
the two approaches. Indeed, the vl28-3 run has a larger amount 
of turbulent energy (Sec l4.2H nside the simulated galaxy clus- 



12 3 4 
AMR eve 

Fig. 4. Ratio of the volume covered by cells at the various AMR 
levels (normalized to the volume of the AMR region), for runs 
V128-3 (red-solid) and dl28-3 (black-dashed) at z = 0.1. The 
thick lines gives the cumulative distribution, while the thin lines 
show the differential distribution. The root grid level is labeled 
as "0". 

ter and this reduces the innermost gas density compared with 
that of d 128; this balances the larger number of cells triggered 
according to their velocity jump. 

4.2. Turbulent motions 

In order to characterize turbulent velocity fields in the com- 
plex environment of galaxy clusters, it is necessary to extract 
velocity fluctuations from a complex distribution of velocities. 
Dolag et al.(2005) proposed that the turbulent gas velocity field 
can be extracted by removing a "local" mean velocity field, 
whose value is obtained by interpolating the 3-D gas velocity 
on large enough scales. Using this approach, it was shown that 
the bulk of laminar infall motions driven by accreted substruc- 
tures in Smoothed Particles Hydrodynamics(SPH) simulations 
develops at scales of the order of > 100 - 300kpc, which in- 
deed corresponds to the core radii of matter clumps accreted by 
massive galaxy clusters. 

Following a similar approach, here we use the ENZO im- 
plementation of the PPM scheme (based on parabolic interpola- 
tions on cells) to map the 3-D local mean velocity field, Vl, and 
for each cell we measure the turbulent velocity as Av = v - Vl; 
v is the gas velocity at the maximum AMR level, while Vl is 
measured at a coarser resolution (for the v256-3 and v256-4 
runs this is A = 292kpc, while for the other runs we consider 
the AMR level corresponding to this scale). We notice that this 
procedure implies a largest possible scale of » 300kpc for tur- 
bulent motions, and therefore in presence of significant turbu- 
lent motions on larger scales our procedure would lead to a 
lower estimate on the total turbulent energy budget. Yet, the 
influence of our filtering scale in the final estimate of the tur- 
bulent energy cannot be larger than a factor ~ 2. This simply 
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comes from the comparison of the kinetic energy and turbu- 
lent energy profiles reported in the last panel of Fig|6] and it is 
consistent with tests previously reported in Dolag et al.(2005) 
e Vazza et al.(2006). The visual inspection (e.g. FigO further 
confirms that most of the velocity structure present in the IGM 
at scales > 300kpc is mostly due to laminar infall motions. 

In all runs, the total mass of the cluster at the center of the 
AMR region at z = is M ~ 2. 1 • 10 14 M©, which corresponds to 
a virial radius of /? V! > = \AMpc. Panels in FigOshow the total 
and turbulent velocity fields for an epoch just after the major 
merger event, z = 0.6, for a slice crossing the AMR region. The 
laminar infall patterns, due to accretion of sub-clumps from 
filaments (left panel), are almost completely removed by our 
filtering of the velocity field, and small scale curling motions 
injected by accreted clumps and by shocks (see also Sec l4.4t 
are well highlighted (right panel). 

The uppermost panels in Fig|6]show the gas density profile and 
the gas entropy profiles of the cluster in all runs. The lower 
panels in the same Figure show the profiles of thermal, tur- 
bulent and kinetic energy, and the ratio between turbulent (or 
kinetic) energy and the total energy E m (kinetic plus thermal) 
inside a given radius. The turbulent energy, E tur b, is measured 
as pAv 2 /2, the total kinetic energy is = pv 2 /2 and the 
thermal energy in the cell is £,/, = (3/2)kspT /pm p ; the ve- 
locity field is always corrected for the galaxy cluster center 
of mass velocity; all profiles refer to z = 0.1. The standard 
AMR run (i.e. over-density based refinement, dl28) shows the 
highest central density and the steepest entropy profile, while 




L : 5 £ 2.5 

\ 1 c^l^ch 1 j 



= 0.6. Right: map of 



all runs with velocity/over-density refinement have flatter pro- 
files. This is explained because merger shocks in runs with the 
velocity/over-density AMR criterion are simulated with higher 
accuracy during cluster lifetime, and they can propagate more 
deeply towards the inner regions of the cluster without be- 
ing damped by resolution effect. At all radii, the runs with 
the velocity/over-density refinement show larger energy bud- 
get in turbulent motions, with a E tur \,jE tot ~ 3 - 4 percent at 
r = 0.1R V i r (E m rb/Eth ~ 5 per cent within the same radius) 
and E tur b/E tot ~ 5 - 8 percent inside R vir (E tur b/E,i, ~ 10-20 
per cent within the same radius). As expected the adoption of a 
larger threshold for 6 (vl28-10) decreases the budget of turbu- 
lent motions in the simulated volume, gradually approaching 
the results of standard AMR(dl28), except for the outermost 
regions, where strong shocks occur and the threshold 6 = 10 
still triggers refinement. Decreasing 6 (vl28-l) increases the 
turbulent energy budget, yet convergence is already reached at 
> 0.2R V i r for 5-3 (vl28-3). The adopting the mesh refinement 
criterion based on velocity jumps for z < 2 (vl28-z2) produces 
profiles consistent with those from the run where this criterion 
is applied since the beginning of the simulation (vl28), pro- 
vided a small difference is found for r ~ 7? V! >. 

In the cases where the AMR peak resolution is fixed at 
A = 36kpc (v256-3,vl28-3,v64-3), the adoption of a larger 
mass resolution in DM particles causes a significant decrease in 
the turbulent budget at large radii (the kinetic energy profiles, 
however, are almost unaffected by that). We find that the rea- 
son for this is that in the cluster outskirts, where strong accre- 



Fig. 5. Left: modulus of total gas velocity in a slice of side 1.5Mpc and depth ISkpc, for the v256-4 run at z 
Mach number (in colors) and turbulent gas velocity field (arrows). 
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Fig. 6. Gas density profiles (top left), gas entropy profiles (top right), E t urb/E t0 ,and Eia n /E tot profiles (lower left) and E, ur b,Eki n 
and E,i ierm profiles (lower right) for all runs of the paper at z = 0.1. The color coding is all the same as in the first panel; E tur b 
refers to gas motions on scales < 3Q0kpc. 



tion shocks are located, satellites with a too coarse DM mass 
resolution have a typically smaller gas and DM density con- 
centration, and they are more easily stripped and/or destroyed 
generating more small scale chaotic motions in the peripheral 
regions of clusters (see also Sec 14.4b . 

The total kinetic energy within 7? V! > in our simulations 
is in line with SPH results with reduced artificial viscosity 
(Vazza et al.2006) and other AMR results obtained with ENZO 
(Iapichino & Niemeyer 2008). However, it is unclear if the 
observed inner turbulence profile can be reconciled with with 
SPH findings, where an increase of the ratio between turbulent 
energy and the total one is observed with decreasing radius, for 
r/R V j r < 0.1 (Dolag et al.2005). On one hand it seems that the 
progressive increase of the DM mass and force resolution in 
our simulations causes the same kind of steepening also in our 
innermost profile, on the other hand the turbulent energy bud- 
get remains smaller by a factor ~ 5 - 6 respect to SPH results. 
Whether or not this is related to the different clusters under 
observation (and to their dynamical states) or if this is this a 
more fundamental issue caused by differences between AMR 
and SPH simulations, is a topic that deserves more accurate 
investigations in the future. 



4.3. Power Spectra and Structure Functions of the 
Turbulent Velocity Field 

We characterize the cluster velocity field through it 3D power 
spectrum, E(k), defined as: 

£(k) = i|v(k)| 2 , (1) 

where v(k) is the Fourier transform of the velocity field: 
v(k) = -i-3 f v(x) e - 2m ' k V*. (2) 

E(k) is calculated with standard FFT algorithm (e.g. 
Federrath, Klessen & Schmidt 2009, and references therein), 
and with a zero-padding technique to deal with the non- 
periodicity of the considered volume . Differently from SPH 
and standard AMR simulations, the velocity plus density re- 
finement allows to follow the cluster velocity field with high 
spatial resolution in lower density regions, with important con- 
sequences on the capability to describe its spectral properties 
over a wide range of scales. 

The left panel in Fig|7]shows the 3-D power spectra calcu- 
lated for all runs at z = 0.1. E(k) is approximately described 
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Fig. 7. Left: 3D power spectra for the velocity field of the various run at z — 0.1. The spectra are shown up to their Nyquist 
frequency; the purple dashed lines shows the -5/3 slope to guide the eye. Right: longitudinal and transverse third-order structure 
functions for velocity field, v (black, and for the density = weighted velocity field, u = p'^v (red), for a sub volume in the vl28-3 
run. The purple dashed line shows the expected scaling for the Kolmogorov model. 



by a simple power law over more than one order of magni- 
tude in k, with a slope not far from a standard Kolmogorov 
model (E(k) oc AT 5 / 3 ). At large scales (k < 4) a flattening in 
the spectrum is observed in all runs, at a wave number roughly 
corresponding to the virial diameter of the cluster, which likely 
identifies the outer scale of turbulent motions connected with 
accretion processes. We remark that for spatial scales < 32A, 
the slope of the power spectrum is affected by the non-uniform 
numerical dissipation that PPM adopts to increase resolution in 
shocks and contact discontinuities (Porter & Woodward, 1994). 
As in the case of the turbulent energy budget, the vl28-10 run 
falls in between the standard AMR run and all the other runs 
with velocity/over-density refinement, while there is almost no 
difference by adopting 6 = 3 or S = 1 as threshold; we find 
no relevant differences if the velocity refinement criterion is 
adopted at z < 2 (vl28-z2) or at z = 30 (vl28). 

Remarkably due to its larger peak resolution, the v256-4 
shows a regular power law for almost two orders of magnitude, 
thus supporting the picture that the simulated IGM is globally 
turbulent starting from sub-Mpc scales. This is also further 
suggested by the right panel in Fig|7] which shows the third or- 
der velocity structure functions for the vl28-3 run, calculated 
as 

V0=<|v(r + D-v(r)| 5 >. (3) 

Shown are the transverse (v ± 1) and longitudinal (v || 1) 
structure functions extracted from a random sub-sample of 
~ 10 5 cells in the simulated volume. For completeness, struc- 
ture functions are also calculated for the density-weighted ve- 
locity, u = p'^v, which was introduced by Kritsuk et al.(2007) 



to study scaling relations for simulated supersonic turbulence. 
All signals show a peak at ~ Mpc scales, thus implying that 
the maximum outer scale that drives turbulence is of the order 
of R vir . 

4.4. Shock Properties 

Shocks in large scale structures have been investigated in a 
number of semi-analytical (Gabici & Blasi 2003; Berrington & 
Dermer 2003; Keshet et al.2003) and numerical works (Miniati 
et al.2001; Ryu et al 2003; Pfrommer et al.2007; Skillman et 
al.2008, Vazza et al.2008; Molnar et al.2009). Observationally, 
merger shocks have been detected only in a few in a few 
nearby X-ray bright galaxy clusters (Markevitch et al.2005; 
Markevitch 2006; Solovyeva et al.2008), and may be possi- 
bly associated with single or double radio relics discovered 
in a number of galaxy clusters (e.g. Roettgering et al.1997; 
Markevitch et al.2005; Bagchi et al.2006; Bonafede et al.2009; 
Giacintucci et al.2008). 

The application of the AMR approach described in this pa- 
per to galaxy allows to follow with high resolution the onset 
and the evolution of shock waves in the IGM of simulated 
galaxy cluster within ~ 2/?,,„- from the clusters center, and to 
explore the connection between shocks and turbulence in large 
scale structures. 

In Fig|8] we present standard 1-D shock-tube tests for a 
weak (M » 1.5) and for a strong (M » 10) shock, where 
we compare the application of the AMR criterion based on ve- 
locity jumps (with 5 = 3), the application of the density jump 
criterion (with 6 P = 2), and a simulation with constant spa- 
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Fig. 8. Top panel: gas density and gas velocity for a 1-D shock 
tube test, at the position of a M « 1.5 shock wave. Bottom 
panel: gas density and gas velocity for a 1-D shock tube test, 
at the position of a M » 10 shock wave. In both cases results 
are shown for a run with N = 1024 cells with uniform resolu- 
tion (unigrid), and for two AMR runs with root grid N = 128 
and 3 additional level of refinement, triggered by density jump 
(d-AMR) and by velocity jump ( v-AMR). The values of gas den- 
sities in both cases have been rescaled by an arbitrary value to 
avoid overlapping with the other lines. 
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Fig. 9. Top: distribution of thermal energy flux at shocks within 
the AMR region of all simulated runs with the over den- 
sity/velocity AMR criterion, at the grid resolution of A = 36 
kpc. Bottom: radial profile of the energy-flux weighted aver- 
age Mach number at z=0.1, for all runs with the over den- 
sity/velocity AMR criterion. 



tial resolution fixed at the maximum resolution level of the 
AMR runs. In both tests, the run with AMR based on veloc- 
ity jumps well matches the results of the fixed resolution run at 
the position of the traveling shock waves, and as result the cor- 
rect shock jump conditions can be basically measured across 3 
cells, in both cases (e.g. Tasker et al.2009). On the other hand 
the AMR based on over density smears the weak shocks across 
a larger distance, since the small density jumps associated with 
M =s 1 .5 is not large enough to trigger any mesh refinement; 
however in the case of the M » 10 shock the gas compres- 
sion is large enough to trigger two level of refinement also in 
the AMR method based on over density. 



Shocks in the our cosmological simulations are identified 
by means of the procedure presented in Vazza, Brunetti & 
Gheller (2009). The algorithm works in the following steps: 

- we consider candidate shocked cells those with V • v < 
(calculated as 3-dimensional velocity divergence to avoid 
confusion with spurious 1-dimensional compressions that 
may happen in very rarefied environments); 

- since shocks in the simulation are typically spread over a 
few cells, we define the shock center with the position of the 
cell in the shocked region with the minimum divergence; 
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- we scan the three Cartesian axes with a one-dimensional 
procedure measuring the velocity jump, Av x y z , between 3 
cells across the shock center; 

- the Mach number of the shock is obtained by inverting 

3 X ' Ml 

where c s is the sound speed in the pre-shock region (the cell 
with the minimum temperature); 

- we finally reconstruct the 3-D Mach number in shocked 
cells as M - (M 2 + M 2 + M 2 ) 1/2 , that would minimizes 
projection effects in the case of diagonal shocks^. 

In order to have a realistic value of the shock Mach number 
at accretion shocks, we apply the post-processing re-ionization 
scheme in Vazza, Brunetti & Gheller (2009), by increasing the 
gas temperature inside cells according to: 

T min (K)=450{^-f w , (5) 
Po 

(where p is the gas density and p = 10~ 32 gr cm 4 ) which 
is found to mimic the integrated effect of the Haardt & Madau 
(1996) run-time re-ionization scheme adopted in ENZO with 
sufficient accuracy. 

The Right panel in Fig|5] shows a map of reconstructed 
Mach number with this method, for a slice of \%kpc crossing 
the AMR region of run v256-4 (with overlaid streamlines of 
the turbulent velocity field). 

In the case that AMR is forced to increase the spatial res- 
olution also around shocks, as our simulations with the novel 
AMR criterion, the Vazza et al.(2009) shock detecting scheme 
can be straightforwardly applied to simulated data at the high- 
est available AMR level (therefore excluding our dl28 run). 
We thus analyze shock statistics in all runs employing the over 
density/velocity AMR criterion at the resolution of A = 36kpc. 

The thermal energy flux across shocks is customary evalu- 
ated as: 

f th =5 M (M)- P M\]l2, (6) 

where p is the pre-shock density and 6m{M) is a monoton- 
ically increasing function of M (e.g. Ryu et al.2003). 

Fig|9](fo/9 panel) shows the distribution of thermal energy 
flux at shocks within the AMR region. The distribution of 
weaker (mostly internal) shocks peaks at M ~ 1.5. Overall, the 
distributions are very steep and consistent with those reported 
in Vazza, Brunetti & Gheller (2009). Compared to Pfrommer 
et al.(2007), who studied shock energetics with high resolution 
GADGET2 simulations, we find significantly steeper energy 
flux distributions in all our runs, a,/, » -3.5 (with /,/,(M)M oc 
M"' h ) for M < 10, compared to a,/, a -2 within the same range 
of Mach number in Pfrommer et al.(2007). The distributions 
of the various re-simulations show relevant differences only for 
shocks with M > 10, where two clear trends can be found: 

1 We notice that recently Skillman et al.(2008) pointed out that the 
measure of Mach numbers with a temperature-based method using a 
coordinate split algorithm, overproduces the number of shocks in the 
case of complex, oblique flows. 



- for a fixed DM mass resolution, the adoption of 6 — 10 
(vl28-10)leads to a significant reduction of the thermal en- 
ergy flux at strong shocks compared to the other runs; 

- for a given AMR criterion based on over density/velocity 
jumps, increasing the DM mass resolution leads to a sig- 
nificant decrease in the thermal energy flux processed at 
strong shocks and to a progressive steepening of the ther- 
mal energy flux distribution. 

Overall we conclude that, once that the velocity jump AMR 
criterion is adopted, the largest amount of difference in the spa- 
tial and energy distribution of shocks is caused by the DM mass 
resolution. 

In order to highlight the reason for this finding, we show in 
Fig[l0]the maps of projected Dark Matter density, gas temper- 
ature, mean Mach number and thermal energy flux for a slice 
of depth=100 kpc for runs v256-3, vl28-3 and v64-3. The in- 
crease in the number of accreted DM clumps in simulations 
with higher DM mass resolution is found to generate a more 
complex temperature distribution, which follows the pattern of 
matter infall on the cluster. On the other hand when the DM 
mass resolution is coarse, outer shocks are found to be more 
regular in shape, and they are characterized by sharper jumps. 
The decrease of DM mass resolution implies that the cluster be- 
comes more spherically symmetric due to the lack of substruc- 
ture, thus our findings qualitatively support those of Molnar et 
al.(2009), which shows that the importance of pressure jumps 
due to accretion shocks in simulated clusters is reduced by a 
factor 5-10 compared to predictions based on spherical models 
(Kocsis et al.2005). 

Fig^\(bottom panel) shows the radial profile of the energy- 
flux weighted average Mach number for all runs with the over 
density/velocity refinement. All runs produce consistent pro- 
files up to R vir , with < M 1.5. As seen above, differences 
are larger at accretion shocks outside Z? VJ >, and in particular we 
find that as soon as the DM resolution is increased, the mean 
strength of shocks at r ~ 1 -2/? v /r 1S reduced by a factor ~ 2-5. 

4.5. Time Evolution 

We produced a highly time-resolved study of turbulence and 
shocks developing in the AMR region of run v256-3. Left panel 
in FigQT]shows the evolution with cosmic time of kE(k) within 
a sub-volume of 3.5Mpc centered on the cluster center. The 
bulk of turbulence injection starts with the onset of the major 
merger, at z ~ 1, and develops at scales in the range ~ 1 - 
2Mpc. At smaller redshifts, the spectrum gradually approaches 
the shape reported in Fig(7] 

To better explore the connection between shock waves and 
turbulence in the merger event, we show in Fis\TT\( right panel) 
the evolution of the thermal energy flux through shocks for the 
same sub-volume considered in the Left panel. The energy flux 
is calculated with the same procedure as in Sec |4.41 but in this 
case no treatment of re-ionization is considered, and therefore 
accretion shocks are stronger than those measured in Fig|9] be- 
cause of the unrealistically low value of gas temperature out- 
side the galaxy cluster at evolved redshifts. 
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A bump of thermal energy flux at strong merger shocks is 
measured at approximately at the same epoch when the bulk 
of large scale kinetic energy is injected in the IGM. Soon after 
virialization occurs, extremely strong shocks become rarer and 
the shocks energy distribution approaches the distribution of 
Left panel in FigJlO] (provided that the considered volume is 
smaller, and that re-ionization is not modeled here). 

5. Conclusions 

A very simple implementation of a new refinement criterion in 
ENZO simulations allows to follow shocks and turbulent mo- 
tions with unprecedented detail up to large distances from clus- 
ter centers. This refinement criterion is successful in catching 
the bulk of turbulent motions developed in the IGM by clus- 
ter formation processes, allows us to measure velocity power 
spectra across two orders of magnitude in spatial scales, and to 
follow shocks statistics and evolution over time in great detail. 

Compared to the standard grid refinement criterion, we find 
that the extra refinement on velocity jumps causes no signif- 
icant extra expense of memory storage, and that by construc- 
tion it readily allows to use accurate shock detecting scheme 
at the largest available resolution in these simulations. In all 
the analyzed runs, the simulated IGM is found to host turbu- 
lent motions (on scales < 300kpc) accounting for a ~ 5 - 25 
per cent of the gas thermal energy within /?,,,>. Compared to re- 
finement based on over-density only, the new criterion shows 
lower inner gas density, flatter entropy profiles, significantly 
larger turbulence budget at all radii and a larger thermal energy 
budget processed at accretion shocks. This is due to the sharper 
representation of shock waves and turbulent motions, and high- 
lights the importance of highly resolving these phenomena in 
discussing accretion processes in the IGM of galaxy clusters. 

When the new over density/velocity AMR criterion is em- 
ployed, the DM mass resolution is found to play a fundamental 
role in setting the properties of the turbulence generation and 
of thermal energy flux at shocks; if DM resolution is increased, 
in-falling matter clumps are less easily destroyed during accre- 
tion and they thus inject less turbulence via the ram pressure 
stripping mechanism. In addition, the complex accretion pat- 
tern established in simulations with high DM mass resolution 
is found to significantly prevent the formation of sharp accre- 
tion shocks, compared to runs where the DM mass resolution 
is coarser. In our simulations we find no relevant differences in 
the properties of turbulence and shock waves if the extra refine- 
ment based on velocity jumps is considered only starting from 
Z < 2. 

Overall, the above results confirm that shocks, turbulence 
and dark matter clustering are inter-playing key ingredients 
which modern cosmological numerical simulations need to fol- 
low with high order accuracy and high resolution to model the 
thermal (and non thermal) properties of the IGM in a realistic 
way. 
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color coding for the liens is shown in the color bar. 
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